Simulation of reconstructions of the polar ZnO (0001) surfaces 
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Surface reconstructions on the polar ZnO (0001) surface are investigated using empirical potential 
models. Several possible reconstructions based around triangular motifs are investigated. The 
quenching of the dipole moment in the material dominates the energetics of the surface patterns so 
that no one particular size of surface triangular island or pit is strongly favoured. We employ Monte 
Carlo simulations to explore which patterns emerge from a high temperature quench and during 
deposition of additional ZnO monolayers. The simulations show that a range of triangular islands 
and pits evolve in competition with one another. The surface patterns we discover are qualitatively 
similar to those observed experimentally. 



I. INTRODUCTION 

Zinc oxide has a wide range of applicability from elec- 
tronics to catalysis *i For example it is used as part of 
ZnO/Me/ZnO (Me=metal) multilayer functional glass 
coatings designed to filter heat-generating infra-red so- 
lar radiation. This is usually achieved by incorporat- 
ing a thin low emissivity metal layer a few nanometres 
thick, 2 In the case where silver is used to construct a 
ZnO-Ag-ZnO sandwich, it has been shown^ that the 
lower Ag(lll)/ZnO(0001) interface may fail for reasons 
not yet fully understood, leading to a sizable cost increase 
in the manufacturing of these devices. To understand this 
effect, one needs to characterise the interfacial structure 
that arises from the growth of Au on Zn(0001). A pre- 
requisite for this is a fundamental understanding of the 
Zn(0001) surface that templates this growth. 

Zinc oxide (zincite) has the well-known wurtzite struc- 
ture with lattice parameters at room temperature and 
ambient pressure of a = 3.25 A, c = 5.207 A, and 
u = 0.3825, and space-group P6smc (no. 186 in crys- 
tallographic tables) The structure may be understood 
as two interpenetrating hexagonal lattices, with each Zn 
(resp. O) sitting at the centre of a distorted O (resp. Zn) 
tetrahedron. The crystal when cut along the (0001) or 
(0001) planes is known to be a type III polar material 
according to the Tasker classification^— That is to say 
that the unit cell is comprised of alternative negative and 
positive charged layers. This ultimately leads to a diverg- 
ing electrostatic potential and should make the two polar 
surfaces of ZnO energetically unfavourable. This, how- 
ever, is not the case as both the O-terminated and Zn- 
terminated polar surfaces show remarkable stability. 10 

Consider a slab of the material with bulk-terminated 
polar surface as used in typical computations (see Figure 
1). Since the polar ZnO(0001) and ZnO(OOl) surfaces 
occur naturally, there must be a mechanism to quench 
the dipole moment that exists normal to the slab surface. 
In order to quench this macroscopic dipole moment a 
transfer of charge across the slab of (1 — 2u)a ~ 0.235 x a 
is necessary, where a is the surface charge density This 
may be understood in terms of the electrostatic energy 
change when charge is moved from one surface to the 




FIG. 1: Bulk crystal structure of wurtzite zinc-oxide with bulk 
lattice parameters a = 3.25 A, c = 5.207 A, and u = 0.3825 A. 
In the bulk each ion is four- fold coordinated, while the surface 
atoms have only three-fold coordination, (colour online) 



other in the direction of the internal electric field of the 
unreconstructed slab. Once sufficient charge has been 
moved, the counter electric field thus established cancels 
the one due to the bulk structure. 

There are several mechanisms which may compensate 
the charge at the surface and counteract the macroscopic 
dipole moment of the semi-infinite crystal. Three, not 
necessarily incompatible, mechanisms have been consid- 
ered in the literature: (i) adsorption of charged species 
e.g. hydroxilation, (ii) modification of the surface region 
by reconstruction, and (Hi) direct charge transfer. 

Until recently, the exact nature of such a mechanism 
in ZnO was not well understood, but recent theoretical 
and experimental studies^— may have resolved the is- 
sue. A combination of surface microscopy techniques and 
Density Functional Theory (DFT) have indicated that, 
depending on the atmospheric environment, mechanisms 
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(i) and (ii) may be active in quenching the polarisation 
of the ZnO(OOOl) surfaces. For the Zn terminated ori- 
entation and in hydrogen-rich conditions, the surface is 
best passivated by adsorption of hydroxyl groups, while 
under low hydrogen partial pressure the surface tends to 
form triangular reconstructions that appropriately com- 
pensate the charge imbalance created by the surface cut. 
This work seems to rule out the third mechanism (hi) 
which had previously been proposed ^involving charge 
transfer between O-terminated and Zn-terminated sur- 
faces. 

The above theoretical studies have relied on DFT 
which, while accurate, has a high computational cost and 
does not allow for a comprehensive search of the phase 
space. For example, in the case of the triangular re- 
construction, STM scans show a range of triangle sizes 
while DFT studies only allow comparison of the ener- 
getics of single, relatively small configurations. In this 
work we present a study which combines fast, albeit less 
accurate, empirical potentials with Monte-Carlo simula- 
tions to study the structure and energetics of the polar 
ZnO(OOOl) surface, focusing on the surface reconstruc- 
tion mechanism to quench the dipole. This approach is 
justified by the dominant role electrostatics plays in the 
surface resconstructionsJ^ 

The rest of the paper has the following structure. Sec- 
tion II describes the methodology employed, discussing 
the empirical potentials, surface relaxation calculations 
and Monte Carlo (MC) simulations. The results are pre- 
sented in Section III, firstly for the energetics of various 
surface reconstructions, and then for the surface patterns 
that emerge from the MC simulations. The implications 
of the results are discussed in the following section, and 
our conclusions are given in the final Section V. 

II. METHODOLOGY 

A. Empirical potentials and surface slab 
calculations 

Previous work by Cat low and coworker s^i^ has shown 
that empirical potentials are well suited to describe the 
details of the structure of the polar surfaces of most ox- 
ides. The parameters for a Buckingham potential are 
fitted following refJ^ to reproduce a variety of properties 
of Zinc Oxide. The potential was kept as simple as possi- 
ble as is appropriate for the desired level of computation. 
For details of the validation of the interatomic potential 
parameters we refer the reader to ref. 17 . The total energy 
is computed by summing all pair interactions of the form 

E(nj) = + Aexp(-r ii /p) - C/r% (1) 

where r^- = ||r^ — Tj\\ is the distance between two ions 
with charges qi and qj . In this work we use formal ionic 
charges ±2e in all our computations. The first term of 
Eq. [T]is the long-range Coulomb pair interaction, while 




FIG. 2: Surface reconstruction \^48 x a/48 with a n = 7 
triangular pit and an additional m — 3 inner pit within the 
larger triangle. The pits and terraces are created by removing 
Zn and O atoms. The topmost layer atoms are shown as large 
spheres (light red for O, dark blue for Zn) and the next layer 
atoms are shown by smaller spheres. The "bulk" atoms in 
the slab (everything below layer 2) are shown only by their 
bonds. Upper left we show the smaller triangular pits in a 
3x3 and 6x6 surface unit cell. 

the second and third terms correspond respectively to the 
repulsive and attractive terms of the short-range Buck- 
ingham pair potential. 

The polarisability effects are described by a core-shell 
model, where the oxygen ion and its electronic cloud are 
modelled by a massive core and a mass-less shell carrying 
different charges (but with total charge — 2e) and linked 
by a spring with energy 

^spring = l_ h5l (2) 

where ki is the spring constant for ion i and 5i is the 
core-shell distance. The empirical parameters A, p, C 
and ki are determined by fitting to available experimen- 
tal properties, such as the elastic constants. The fitting 
was performed with the GULP code 18 using 8 potential 
parameters and various parameters from observable data. 
The computed bulk properties are compared to some rel- 
evant experimental values in Table HI and the values of 
the parameters used in this work are given in Table [TTJ 

Surface structure calculations were performed using 
three-dimensional periodic slabs with a large vacuum 
gap normal to the (0001) surface. Each slab contains 
one Zn-terminated surface and one O-terminated surface 
(see Fig. 1). Several reconstructions were created at the 
surfaces, where overall charge neutrality was ensured by 
removing the oppositely charged species from the other 
side of the slab. The outermost three surface layers on 
each side of the slab were allowed to relax. The energy 
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TABLE I: Comparison of some of the computed bulk prop- 
erties to available experimental data. The experimental data 
are taken from standard tables 19 . 



ZnO (wurtzite) This work Experiment 



a(A) 


3.27 


3.25 


c(A) 


5.18 


5.207 


U 


0.3819 


0.3825 




4.22 


9.26 


£33 


4.59 


11.0 


Cn(GPa) 


222.22 


209.7 


C 33 (GPa) 


220.14 


210.9 



of different surface reconstructions were computed and 
an energy hierarchy constructed by comparing their sur- 
face energies. In order to estimate the surface energy 
of a given surface structure we use the total energy of 
the bulk unit cell as a reference state, in which case the 
surface energy is given by 

7 = (Esiab ~ N ce u x e hulk ) /2A (3) 

where A = (N s a) 2 sin 7 is the surface area of a slab (see 
Fig. [T]) containing N s x N s x N z lattice units, £^ a 6 is 
its relaxed total energy, and £buik is the energy of a bulk 
unit cell. The quantity N ce u is the effective number of 
unit cells in the slab calculated by dividing the number of 
atoms in the slab by the number of atoms in a bulk unit 
cell ( four in the case of ZnO). The factor of two accounts 
for the fact that the surface energy is the average surface 
energy of both sides of the slab. 

It is important to emphasise again that the surface en- 
ergy calculated in this way depends on slab thickness, 
unless we have perfect quenching of the dipole moment 
by creating a net charge transfer of 0.235 x a from one sur- 
face to the other in the reconstruction. For this reason, 
we will compare the energies of various surface structures 
using the same slab thickness N z = 6 (twelve bilayers). 

B. Monte-Carlo simulations 

While a large number of structures may be explored 
using the above empirical model, it is impractical to find 
the lowest energy reconstruction using more and more 
elaborate guesses of the surface structure. In order to 
explore the large phase space of possible surface recon- 
structions we used Monte Carlo (MC) simulations with 
bulk lattice positions in a slab. The ions in the three 
upper-most bilayers of the slab are allowed to hop within 
their own bilayer and into the bilayers directly above 
or below. The simulations were started from different 
initial configurations and the system left to evolve ac- 
cording to the Metropolis algorithm. The initial con- 
figurations of the three uppermost bilayers were formed 
either by ion removals or by addition of ions to the clean 
slab. The overall charge neutrality was again ensured by 
adding/removing the oppositely charged species from the 



TABLE II: Interatomic potential parameters used for the 
Buckingham potential. These were obtained by fitting to the 
experimental parameters of Tab. The spring Constants are in 
eV.A- 2 : ko = 15.52 , k Zn = 8.57. 

A(eV) p(A) C(eVxA b ) r cutoff 
Zn-O s 499.6 0.359 O0 0-10 A 
O s -O s 22764.0 0.149 27.88 0-12 A 

other side of the slab at the start of each simulation. As 
several studies^ have shown that the surface layer relax- 
ation is less than 0.1 A, for the sake of simplicity (and 
computational efficiency) we have neglected the effect of 
lattice relaxation on surface energy during these MC sim- 
ulations. 

A periodic slab model is used throughout, with N z = 6 
bulk unit cells along the c axis and a large vacuum of 
L z = 30A added to form the super-cell. We verify that 
L z is large enough by ensuring that the surface energy 
does not depend on L z . 

Successive configurations are generated by a series of 
nearest neighbour hops of either species selected at ran- 
dom at the uppermost Zn-terminated surface, with ac- 
ceptance probability e - AE / k sT w j iere j 1 \ s th e temper- 
ature and AE is the difference between total energies 
of the successive trial configurations. The bottom O- 
terminated surface reconstruction remains fixed in the 
simulations. The energy of a given atom in the slab is 
simply the sum of all its pair interaction with the other 
ions in the slab, with the pair interaction of ion i given 

by 

Si = Y^E( rij ). (4) 

3 

Here, also for computational efficiency, we neglect the 
shell model component of the potential E spring in the 
MC work only. 

Since in this system we have no mechanism to quench 
any dipole across the slab, we ensure that the simulation 
is started from configurations with only a small residual 
dipole. The initial configuration of the MC simulation is 
then disordered by running the simulation at very high 
temperature leading to a fully disordered arrangement of 
the surface species, after which the temperature is low- 
ered abruptly. The simulations are run at high temper- 
ature for a large enough number of steps that the initial 
ordering disappears. 

The most expensive step of the simulation is the en- 
ergy evaluation which includes long-range terms. The 
Coulomb sum being conditionally convergent in a peri- 
odic system, we make use of the Ewald sum2£ — 

l rv erfc(ar^ n ) 
B « = oL — (5) 

* . . I ij n 

^ = yEV |s(k)|2 (6) 
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TABLE III: Summary of the surface energies of several slabs, 
all with N z = 6, and various total number of ions N. For 
the relaxed structures, the total energies are computed using 
GULP with the potential parameters listed in table [III The 
three outermost bilayers are allowed to relax while the rest of 
the slab is kept fixed. The vacancy concentration Gvac refers 
to the ratio of excess zinc vacancies in the uppermost surface 
layer created by removing zincs and oxygen from both sides 
of the slab. 



Ovac(ML) 



~N -y{eV/A*) 
rigid relaxed 



Bulk 

2x2 
4x4 
I x a/48 
6x6 



E 



0.500 
0.250 
0.208 
0.055 



88 4.210 
352 0.255 
1034 0.306 
856 2.248 



4.162 
0.086 
0.245 
1.447 



2tt 

V 



M 



(7) 



where qi is the formal charge of on z, r^is the position 
of the ion within the periodic slab, S(k) = qie lVxi is 
the structure factor, and M z is the z coordinate of the 
total dipole moment in the slab M = The pa- 

rameter a is determined using the requirement that the 
Ewald sum is accurate yet efficient (see for example^). 
It is worth noting that this expression is the more com- 
putationally efficient 3D version of the Ewald sum, not 
the two-dimensional version. If the vacuum slab is cho- 
sen large enough, only a correction due to the residual 
surface dipole is necessary. 

The above total energy is computed once at the be- 
ginning of the run, and updated in the course of the 
simulation by only computing the energy difference be- 
tween trial configurations. This considerably speeds up 
the computation of the energy and scales as TV 1 / 2 , where 
N is the number of particles in the system. 



III. RESULTS 

A. Surface reconstructions and energy hierarchy 

On the zinc-terminated surface, it was experimentally 
shown that the triangular reconstructions are a single 
bilayer high. We therefore form the reconstructions by 
removing Zn and O atoms from the top-most bilayer and 
refer in what follows to the excess zinc vacancy concen- 
tration Q vac - Recall that we remove the opposite charge 
species from the O- terminated surface, keeping overall 
charge neutrality and imposing an electric field across 
the slab which compensates the bulk dipole moment for 
©vac = 0.235. The stability of the surface reconstruc- 
tion is assessed by studying the surface energy as a func- 
tion of this excess Zn vacancy concentration. Table IIIII 
shows the surface energies for various zinc vacancy con- 
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FIG. 3: Change in surface energy as a function of excess zinc 
vacancy on the surface expressed in fractions of a monolayer, 
where A7 = ^fnvac — Jdean is the net surface energy compared 
to the surface energy of the bulk-terminated surface. The blue 
line refers to bulk position of the atoms while the red curve 
correspond to optimised position for the three bilayers nearest 
the surface. The surface energy change when creating isolated 
vacancies is shown by the black curve. Note that the small 
spheres refer to the topmost bilayer, while the larger ones 
refer to the bilayer immediately below it. 



centrations and their corresponding surface reconstruc- 
tions illustrated in Figure [2] . Note that the reconstruc- 
tion in the \/48 x \/48 surface unit cell as shown in Fig. 
[2] is not the most stable in this work, contrary to the 
prediction from density functional theory. This may be 
explained by a simple electrostatic argument, where the 
formal ionic representation of the species tends to overes- 
timate the contribution of the Coulomb interaction, and 
thus favours smaller surface reconstructions. 

We have computed the surface energies of several more 
surface reconstructions, for varying values of the excess 
zinc concentration, as well as for isolated vacancies. The 
results are summarised in Fig. [3] and are plotted in ref- 
erence to the surface energy of the bulk-terminated slab 
with N z = 6. We reiterate that the surface energy of the 
bulk-terminated slab is ill-defined due to the presence of 
a large surface dipole. It is only the difference in surface 
energies of various reconstructions that is of interest. The 
most stable reconstruction is the small 4x4 cluster with 
the smallest triangular vacancy (see Figure [2]), with the 
more exotic reconstruction on the \/48 x \/48 surface cell 
lying nearby. The isolated vacancies consistently have a 
larger surface energy. Relaxation of the surface layers 
has a small effect on the surface stability, but does not 
impact the overall ordering of the various reconstruction 
energies. 
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In conclusion, the surface hierarchy obtained using the 
simple empirical potential shows no preference for large 
triangular reconstructions. Instead, small clusters which 
locally quench the surface dipole are preferred. The small 
impact of surface relaxation reflects the dominance of 
dipole moment quenching. This means that surface re- 
laxation can be neglected in the MC simulations which 
follow, in which we allow the system to explore configu- 
ration space to see if larger reconstructions emerge nat- 
urally in our model with formal ionic charges. 



B. Monte-Carlo Simulations 

We have performed the MC simulation of a high tem- 
perature quench on a slab with N z = 6, N s =32 and 
with the surfaces initially tessellated with the 4x4 tri- 
angular reconstructions giving N = 22, 528 ions in total. 
Representative snapshots of the species distribution in 
the three uppermost layers of the Zn-terminated surface 
are shown in Fig. H^a-c). We show the three uppermost 
bilayers separately (bilayer 1-3) at (a) the start of the 
simulation, (b) at the end of the high temperature run, 
and (c) after quenching. The figure clearly shows the 
effect of the high temperature, with the ordered trian- 
gular reconstructions completely disappearing from the 
first bilayer. The second and third bilayers are occupied 
as well, with no apparent ordering. After the system is 
cooled and left to evolve, large triangular reconstructions 
in the lowest bilayer begin to nucleate via aggregation of 
small triangular units. The smallest such unit is the one 
unit-cell triangle (three O ions surrounding one Zn ion) 
and is consistent with the earlier empirical potential re- 
sult. The second bilayer only shows a few isolated a— side 
triangles, while the third layer is now completely empty. 
These observations are qualitatively identical for a large 
class of system size and parameters. 

In Fig. 5 we show the evolution of the surface energy 
and the dipole moment normal to the surface in the sim- 
ulation. During the high temperature anneal, both the 
surface energy and dipole moment are large in magni- 
tude. This is due to the almost random placement of the 
surface layer ions into the 3 accessible bilayers, which 
leads to an obvious loss of bonding energy. The dipole 
moment also increases in magnitude since there is an ex- 
cess of oxygen ions over zinc in the surface, so displacing 
them on average by one bilayer changes the total moment 
in the system. Upon quenching, the surface energy and 
dipole moment quickly reduce in magnitude, and in fact 
reach slightly lower values than in the starting configura- 
tion. At the start of the simulation, the dipole moment 
is —162.29 e.A, and is not zero since O vac = 0.25 rather 
than the ideal & vac = 0.235. After the temperature is 
raised, the dipole moment is roug hly -300 e.A, since the 
topmost ions are now randomly distributed. By the end 
of the quench, after Nmc = 4 x 10 5 , its value is slightly 
lower, —236.15 e.A. This is achieved by the few ions oc- 
cupying the second bilayer. For the surface energy, the 



quenched value is 0.36 eV/A 2 , as opposed to 0.59 eV/A 2 
after heating. 

The quenched MC simulation clearly shows that the 
system can evolve to a structure with preferred triangu- 
lar motif and no regular tessellation of the surface. It is 
of interest to see if this effect also emerges during a simu- 
lation where the total number of ions in the surface layers 
increases over time, mimicking epitaxial growth. Starting 
from the end of a quenched simulation, we perform the 
MC annealling but now add pairs of Zn 2+ /0 2_ ions at 
separate, randomly chosen sites in the three uppermost 
layers of the slab. Results from a simulation performed 
on a N z = 6, N s = 32 slab (N = 8, 800 ) are shown 
in Fig. [6] Again, the simulation shows that triangu- 
lar reconstructions form spontaneously, and grow larger 
by nucleation from characteristic smaller aggregates. We 
also see that the growth on the 2nd bilayer proceeds be- 
fore the 1st bilayer is complete, leading to a surface with 
multiple ad-islands and pits of various size. 

In Fig. [5] we also show how the surface energy and 
dipole moment change during the deposition simulation. 
The addition of the ions allows the system to find struc- 
tures with decreasing magnitude of dipole moment, since 
thereby lowering the electrostatic energy. An interest- 
ing feature shown by Fig. [5] is that both surface energy 
and dipole moment oscillate as more ions are added, and 
reach even lower values as we deposit more and more ions. 
Thus the system finds a steady state with lower energy 
and dipole moment by making use of the larger number 
of degrees of freedom made available by the deposited 
ions. 



IV. DISCUSSION AND CONCLUSION 

It is observed from the STM results^ that not one 
but several triangular features co-exist at the surface of 
Zn-terminated ZnO(0001). This behaviour also emerges 
in our Monte-Carlo simulations, for both the quenched 
structure in Fig. 2] and the deposition structure in Fig. 
[6j The reason for this can be traced to the small energy 
differences between the various surface reconstructions 
shown in Table 1111} there is no substantially preferred 
reconstruction, provided # V ac is close to 0.235 locally. 
Therefore the patterns that emerge in the MC simula- 
tions result from the competitive growth of energetically 
comparable triangular reconstructions. For this reason, 
we do not observe any long time coarsening of the struc- 
tures in the simulation, as can be seen from the surface 
energy evolution shown in Fig. [5] during the quenched 
phase before deposition starts, and neither do we ob- 
serve the system being restored to its regular tessellated 
starting configuration. 

Our quenched MC simulation in Fig. H] produces an 
interesting surface morphology that is reminiscent of the 
experimental STM image g 11 ! 13 . Furthermore, the sim- 
ulation with increasing surface coverage of Fig. [6] also 
includes other morphological features such as co-existing 




FIG. 4: Snapshot of Zn 2+ and 2 " ions in the three topmost layers at different stages of the evolution of the MC simulation. 
The initial configuration (a) is evolved at high temperature until the initial reconstructions are melted giving rise to a random 
distribution of species in all three layers (b). The temperature is then abruptly at which point triangular structures begin to 
spontaneously in bilayer 1, while the top bilayers are gradually emptied, (c). 



ad-islands and pits that are found experimentally. There- 
fore, we believe that the simulations capture some of the 
main physical processes that give rise to these surface 
reconstructions. However, it is important to note that 
the reconstructions found by the simulation are much 
smaller than those observed experimentally. The largest 
triangle observed from the MC simulation has a side of 
the order of 20Athus only reproducing the smallest clus- 
ters observed in the STM scans. It is possible that this 
is due to the lattice sizes and simulation durations that 
are accessible. Another explanation is that our use of 
formal charges (=L2e) on the ions is not justified in this 
system. Using formal charges probably overestimates the 
strength of the Coulomb interaction, thereby tending to 
make the triangular clusters at the surface more com- 
pact. Work to augment the current model using a charge 



equilibration (or QEq) scheme^ is therefore planned. 

In conclusion, we believe that the models presented 
here do help to explain the Zn-terminated surface recon- 
structions observed experimentally. Whilst the accuracy 
of our models cannot compete with DFT, we are not 
restricted to studying individual reconstructions. The 
use of empirical potentials allows us to explore the phe- 
nomenology of the surface reconstructions more freely, 
and we find that a broad range of characteristic triangu- 
lar motifs naturally emerge in our simulations, qualita- 
tively consistent with the STM results. 
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FIG. 6: In a way similar toHJwe show snapshots of the three upper bilayers as the deposition of a single monolayer (ML) is 
carried out. The top row shows the three bilayers before the deposition is started, the middle one after 0.5 ML was deposited, 
and the lower one after a full ML was deposited onto the surface. Note that the triangular features of the reconstruction are 
essentially conserved throughout the deposition. 



